Analysis parameters:
iGraph Parameters
Load intermediate data after running
worksapce_ml.RData
load(file = "../data/__cache__/worksapce_ml.RData")
USE_NESTED <- TRUE
SAVE_PLOT <- FALSE
To calculate the averaged corr matrix A
C.z <- FisherZ(C)
C.zmean <- matrix(colMeans(C.z), nrow=264, ncol = 264)
A <- FisherZInv(C.zmean)
A.vec <- as.vector(A)
Calculate W from betas and A Matrix
connectom2matrix <- function(connectome, w) {
empty_mat <- matrix(0, 264, 264, dimnames = list(paste0("X", 1:264), paste0("X", 1:264)))
empty_mat[connectome$index] <- connectome$W
return(empty_mat)
}
# convert a 264*264 matrix back to connectom df
matrix2connectom <- function(mat, connectome, col_name) {
connectome$temp = mat[connectome$index]
connectome <- rename(connectome, !!col_name := temp)
return(connectome)
}
# make the matrix symmetric
make_symmetric <- function(m) {
# lower.tri is 0.0
m[lower.tri(m)] <- t(m)[lower.tri(m)]
return(m)
}
power_atals <- power2011 %>%
rename(ROI.Name = ROI, x.mni=X, y.mni=Y, z.mni=Z, network=NetworkName) %>%
mutate(ROI.Name=as.integer(ROI.Name), index = as.integer(ROI.Name),
x.mni=as.integer(x.mni), y.mni=as.integer(y.mni), z.mni=as.numeric(z.mni)) %>%
dplyr::select(ROI.Name, x.mni, y.mni, z.mni, network, index)
check_atlas(power_atals)
if(USE_NESTED) {
connectome_data <- connectome_nested
} else {
connectome_data <- connectome
}
Wconnectome <- connectome_data %>%
mutate(A = A.vec[connectome_data$index], W = A*Beta)
#separate(connection, c("connection1", "connection2"))%>%
#separate(network, sep = "-", c("network1", "network2"), remove = F) %>%
#filter(str_detect(network, pattern = "-1-")) %>%
#network1 = ifelse(str_detect(network, pattern = "-1-"), -1, network1)) %>%
#mutate(connection_type = ifelse(network1==network2, "Within", "Between"))
W_mat <- matrix(0, ncol = 264, nrow = 264)
W_mat[Wconnectome$index] <- Wconnectome$W
W_mat <- make_symmetric(W_mat) #CHECKED correct W_mat
rownames(W_mat) = power_atals$network
colnames(W_mat) = power_atals$network
power_color = power2011 %>%
filter(NetworkName!="Sensory/somatomotor Mouth") %>%
select(NetworkName, Color) %>%
distinct()
colors <- c(Uncertain = power_color$Color[9], #power_color$Color[1],
`Sensory/somatomotor Hand` = power_color$Color[6], #power_color$Color[2],
`Cingulo-opercular Task Control` = power_color$Color[3],
Auditory = power_color$Color[4],
`Default mode` = power_color$Color[5],
`Memory retrieval?` = power_color$Color[2], #power_color$Color[6],
`Ventral attention` = power_color$Color[7],
Visual = power_color$Color[8],
`Fronto-parietal Task Control` = power_color$Color[1], #power_color$Color[9],
Salience = power_color$Color[10],
Subcortical = power_color$Color[11],
Cerebellar = power_color$Color[12],
`Dorsal attention` = power_color$Color[13])
colors <- c(Uncertain = power_color$Color[1],
`Sensory/somatomotor Hand` = power_color$Color[2],
`Cingulo-opercular Task Control` = power_color$Color[3],
Auditory = power_color$Color[4],
`Default mode` = power_color$Color[5],
`Memory retrieval?` = power_color$Color[6],
`Ventral attention` = power_color$Color[7],
Visual = power_color$Color[8],
`Fronto-parietal Task Control` = power_color$Color[9],
Salience = power_color$Color[10],
Subcortical = power_color$Color[11],
Cerebellar = power_color$Color[12],
`Dorsal attention` = power_color$Color[13])
if (SAVE_PLOT) {
png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure1.png", bg = "transparent",
width = 800 , height = 800,
units = "px", res = 150)
circos.clear()
chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1,
symmetric = TRUE, scale = TRUE, reduce = FALSE,
annotationTrackHeight = mm_h(c(15, 1)),
annotationTrack = c("grid"),
grid.col = colors, col = ifelse(W_mat>0, "tomato", "#00000000"))
title("Predictive Group: Declarative (W > 0)", outer = F, cex.main = 1.5)
} else {
circos.clear()
chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1,
symmetric = TRUE, scale = TRUE, reduce = FALSE,
#annotationTrackHeight = mm_h(c(15, 1)),
#annotationTrack = c("grid"),
grid.col = colors, col = ifelse(W_mat>0, "tomato", "#00000000"))
title("Predictive Group: Declarative (W > 0)", outer = F, cex.main = 1.5)
}
if (SAVE_PLOT) {
png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure2.png", bg = "transparent",
width =800 , height = 800,
units = "px", res = 150)
circos.clear()
chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1,
symmetric = TRUE, scale = TRUE, reduce = FALSE,
annotationTrackHeight = mm_h(c(15, 1)),
annotationTrack = c("grid"),
grid.col = colors, col = ifelse(W_mat<0, "steelblue", "#00000000"))
title("Predictive Group: Procedural (W < 0)", outer = F, cex.main = 1.5)
} else {
circos.clear()
chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1,
symmetric = TRUE, scale = TRUE, reduce = FALSE,
#annotationTrackHeight = mm_h(c(15, 1)),
#annotationTrack = c("grid"),
grid.col = colors, col = ifelse(W_mat<0, "steelblue", "#00000000"))
title("Predictive Group: Procedural (W < 0)", outer = F, cex.main = 1.5)
}
circos.clear()
chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1,
symmetric = TRUE, scale = TRUE, reduce = FALSE,
order = unique(power2011$NetworkName),
preAllocateTracks = 1,
annotationTrackHeight = mm_h(c(10, 10)),
annotationTrack = c("grid"),
grid.col = colors, col = ifelse(W_mat>0, "tomato", "steelblue"))
title("Declarative + Procedural Network", cex.main = 2)
circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
xlim = get.cell.meta.data("xlim")
ylim = get.cell.meta.data("ylim")
sector.name = get.cell.meta.data("sector.index")
circos.text(mean(xlim), ylim[1] + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))
circos.axis(h = "top", labels.cex = 0.5, major.tick.percentage = 0.2, sector.index = sector.name, track.index = 2)
}, bg.border = NA)
Next, we will look at Graph properties of two networks
###Graph Density
The proportion of present edges from all possible edges in the network.
# select cols
roi_links <- Wconnectome %>% dplyr::select(connection1, connection2, W, connection_type, network, network_names)
# rename cols
colnames(roi_links) <- c("from", "to", "weight", "connection_type", "network", "network_names")
roi_nodes <- power2011 %>% rename(id = ROI) %>%
mutate(NetworkName=factor(NetworkName),
Color=factor(Color))
levels(roi_nodes$Color) <- sample(colors(T), 14)
# create a graph
net <- graph_from_data_frame(d=roi_links, vertices=roi_nodes, directed=F)
#g <- graph_from_adjacency_matrix(W_mat,, mode = "upper")
net.d <- net - E(net)[E(net)$weight<0]
net.p <- net - E(net)[E(net)$weight>0]
df.density <- data.frame("edge_density"=c(edge_density(net.d, loops=F),
edge_density(net.p, loops=F),
edge_density(net, loops=F)),
"network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))
df.density %>% ggbarplot(x="network", y="edge_density", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
theme_pander() +
ggtitle("Degree Edge Density")
A network diameter is the longest geodesic distance (length of the shortest path between two nodes) in the network. In igraph, diameter() returns the distance, while get_diameter() returns the nodes along the first found path of that distance.
# make negative weights to positive
net.p.abs <- net.p
E(net.p.abs)$weight <- E(net.p.abs)$weight * (-1)
# make negative weights to positive
net.abs <- net
E(net.abs)$weight[E(net.abs)$weight<0] <- E(net.abs)$weight[E(net.abs)$weight<0] * (-1)
df.diameter <- data.frame("diameter"=c(diameter(net.d, directed=F),
diameter(net.p.abs, directed=F),
diameter(net.abs, directed=T)),
"network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))
df.diameter %>% ggbarplot(x="network", y="diameter", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
theme_pander() +
ggtitle("Network Diameter")
Centrality functions (vertex level) and centralization functions (graph level). The centralization functions return res - vertex centrality, centralization, and theoretical_max - maximum centralization score for a graph of that size. The centrality function can run on a subset of nodes (set with the vids parameter). This is helpful for large graphs where calculating all centralities may be a resource-intensive and time-consuming task.
Centrality is a general term that relates to measures of a node’s position in the network. There are many such centrality measures, and it can be a daunting task to wade through all of the different ways to measure a node’s importance in the network. Here, we will introduce just a few examples.
df.centrality <- data.frame("centr_degree"=c(centr_degree(net.d, normalized=T)$centralization,
centr_degree(net.p, normalized=T)$centralization,
centr_degree(net, normalized=T)$centralization),
"network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))
df.centrality %>% ggbarplot(x="network", y="centr_degree", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
theme_pander() +
ggtitle("Network Degree Centrality ")
Let’s now do the same for betweenness centrality, which is defined as the number of geodesic paths (shortest paths) that go through a given node. Nodes with high betweenness might be influential in a network if, for example, they capture the most amount of information flowing through the network because the information tends to flow through them. Here, we use the normalized version of betweenness.
Closeness (centrality based on distance to others in the graph) Inverse of the node’s average geodesic distance to others in the network.
df.sloseness <- data.frame("centr_clo"=c(centr_clo(net.d)$centralization,
centr_clo(net.p)$centralization,
centr_clo(net)$centralization),
"network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))
df.sloseness %>% ggbarplot(x="network", y="centr_clo", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
theme_pander() +
ggtitle("Network Closeness")
df.distance <- data.frame("mean_distance"=c(mean_distance(net.d, directed=F),
mean_distance(net.p, directed=F),
mean_distance(net, directed=F)),
"network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))
df.distance %>% ggbarplot(x="network", y="mean_distance", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
theme_pander() +
ggtitle("Network Closeness")
df.assortativity<- data.frame("assortativity_degree"=c(assortativity_degree(net.d, directed=F),
assortativity_degree(net.p, directed=F),
assortativity_degree(net, directed=F)),
"network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))
df.assortativity %>% ggbarplot(x="network", y="assortativity_degree", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
theme_pander() +
ggtitle("Network Assortativity Degree")
> Combined Table
options(scipen=999)
df.density %>%
select(network, edge_density) %>%
left_join(df.diameter) %>%
left_join(df.centrality) %>%
left_join(df.sloseness) %>%
left_join(df.distance) %>%
left_join(df.assortativity) %>%
rename("centrality_degree"=centr_degree,
"centrality_closeness"=centr_clo) %>%
kable(format.args = list(scientific = TRUE, big.mark = ",", digit=3))
| network | edge_density | diameter | centrality_degree | centrality_closeness | mean_distance | assortativity_degree |
|---|---|---|---|---|---|---|
| Declarative | 9.22e-04 | 1.12e+00 | 6.68e-03 | 5.00e-05 | 1.14e+00 | -1.85e-01 |
| Procedural | 1.56e-03 | 2.73e+00 | 1.37e-02 | 1.83e-04 | 1.83e+00 | -8.19e-02 |
| Full | 2.48e-03 | 2.73e+00 | 1.27e-02 | 3.04e-04 | 2.56e+00 | 4.60e-02 |
#colors <- factor(power_atals$network)
#levels(colors) <- colors14
#power_atals$colors <- as.character(temp)
check_atlas(power_atals)
x1 <- W_mat
x1[x1<0] <- 0
p1 <- brainconn(atlas=power_atals, conmat=x1, node.color = power2011$Color, view = "top",
node.size = igraph::degree(net.d)*2.5, all.nodes = TRUE,
edge.color = "tomato", edge.color.weighted = FALSE, scale.edge.width = c(1,3), edge.alpha = 0.6,
label.edge.weight = F, show.legend = F,
background.alpha = .3) #+ ggtitle("Strategy Predictability: W")
p2 <- brainconn(atlas=power_atals, conmat=x1, node.color = power2011$Color, view = "left",
node.size = igraph::degree(net.d)*2.5, all.nodes = TRUE,
edge.color = "tomato", edge.color.weighted = FALSE, scale.edge.width = c(1,3), edge.alpha = 0.6,
label.edge.weight = F, show.legend = F,
background.alpha = .3) #+ ggtitle("Strategy Predictability: W")
p3 <- brainconn(atlas=power_atals, conmat=x1, node.color = power2011$Color, view = "back",
node.size = igraph::degree(net.d)*2.5, all.nodes = TRUE,
edge.color = "tomato", edge.color.weighted = FALSE, scale.edge.width = c(1,3), edge.alpha = 0.6,
label.edge.weight = F, show.legend = F,
background.alpha = .3) #+ ggtitle("Strategy Predictability: W")
x2 <- W_mat
x2[x2>0] <- 0
p4 <- brainconn(atlas=power_atals, conmat=x2*-10, node.color = power2011$Color, view = "top",
node.size = igraph::degree(net.p)*2.5, all.nodes = TRUE,
edge.color = "steelblue", edge.color.weighted = FALSE,
scale.edge.width = c(1,3), edge.alpha = 0.6,
label.edge.weight = F, show.legend = F,
background.alpha = .3)
p5 <- brainconn(atlas=power_atals, conmat=x2*-10, node.color = power2011$Color, view = "left",
node.size = igraph::degree(net.p)*2.5, all.nodes = TRUE,
edge.color = "steelblue", edge.color.weighted = FALSE,
scale.edge.width = c(1,3), edge.alpha = 0.6,
label.edge.weight = F, show.legend = F,
background.alpha = .3)
p6 <- brainconn(atlas=power_atals, conmat=x2*-10, node.color = power2011$Color, view = "back",
node.size = igraph::degree(net.p)*2.5, all.nodes = TRUE,
edge.color = "steelblue", edge.color.weighted = FALSE,
scale.edge.width = c(1,3), edge.alpha = 0.6,
label.edge.weight = F, show.legend = F,
background.alpha = .3)
if (SAVE_PLOT) {
png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure4.png", bg = "transparent",
width =500, height = 500, units = "px", res = 150) }
p1
if (SAVE_PLOT) {
png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure5.png", bg = "transparent",
width =500, height = 500, units = "px", res = 150)}
p2
if (SAVE_PLOT) {
png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure6.png", bg = "transparent",
width =500, height = 500, units = "px", res = 150)}
p3
if (SAVE_PLOT) {
png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure7.png",
bg = "transparent", width =500, height = 500, units = "px", res = 150)}
p4
if (SAVE_PLOT) {
png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure8.png", bg = "transparent",
width =500, height = 500, units = "px", res = 150)}
p5
if (SAVE_PLOT) {png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure9.png", bg = "transparent",
width =500, height = 500, units = "px", res = 150)}
p6
Look at the distribution of network in two groups
DUPLICATE <- TRUE
c1 <- Wconnectome %>% filter(W>0) %>% mutate(roi = as.integer(connection1)) %>% dplyr::select(roi) %>% unlist()
c2 <- Wconnectome %>% filter(W>0) %>% mutate(roi = as.integer(connection2)) %>% dplyr::select(roi) %>% unlist()
c3 <- Wconnectome %>% filter(W<0) %>% mutate(roi = as.integer(connection1)) %>% dplyr::select(roi) %>% unlist()
c4 <- Wconnectome %>% filter(W<0) %>% mutate(roi = as.integer(connection2)) %>% dplyr::select(roi) %>% unlist()
if (DUPLICATE) {
c12 <- c(c1, c2)
c34 <- c(c3, c4)
} else {
c12 <- unique(c(c1, c2))
c34 <- unique(c(c3, c4))
}
df.c1 <- power2011[c12,] %>%
mutate(NetworkName = factor(NetworkName)) %>%
count(NetworkName, name = "count", .drop = F)%>%
right_join(power2011 %>% dplyr::select(NetworkName, Color) %>% distinct(), on="NetworkName") %>%
mutate(count=as.integer(ifelse(is.na(count), 0, count))) %>%
arrange(NetworkName)
p7 <- power_color %>%
left_join(df.c1) %>%
ggplot(aes(x=count, y=NetworkName)) +
geom_col(fill = power_color$Color) +
scale_fill_manual(guide = guide_legend(reverse = F),
values = df.c1$Color,
labels = df.c1$Color) +
scale_x_reverse() +
theme_pander() +
ggtitle("Distribution of connections", subtitle = "Declarative Network") +
theme(legend.position = "right",
axis.text.y = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size = 20),
axis.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.text = element_text(size = 20))
p7
# count number of networks included in
df.c2 <- power2011[c34,] %>%
mutate(NetworkName = factor(NetworkName)) %>%
count(NetworkName, name = "count", .drop = F) %>%
right_join(power2011 %>% dplyr::select(NetworkName, Color) %>% distinct(), on="NetworkName") %>%
mutate(count=ifelse(is.na(count), 0, count)) %>%
arrange(NetworkName)
p8 <- power_color %>%
left_join(df.c2) %>%
ggplot(aes(y = count, x=NetworkName)) +
geom_col(fill = power_color$Color) +
coord_flip() +
scale_fill_manual(guide = guide_legend(reverse = F), values = power_color$Color) +
theme_pander() +
ggtitle("Distribution of connections", subtitle = "Procedural Network") +
theme(legend.position = "right",
axis.text.y = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size = 20),
axis.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.text = element_text(size = 20))
p8
if (SAVE_PLOT) {
png(filename ="analysis_lasso_re_3connectom/figure-html/network_figure10.png", bg = "transparent",
width = 2000, height = 2000,
units = "px", res = 300)}
ggbarplot(df.c2 , x="NetworkName", y="count", fill = "NetworkName", color="white",
palette = df.c2$Color, #order = "count",
legend = "right",
rotate = TRUE, ggtheme = theme_pander(),
title = "Distribution of connections") +
labs(fill = "Network Name")
# ggbarplot(df.c2 , x="NetworkName", y="count", fill = "NetworkName", color="white",
# palette = df.c2$Color, #order = "count",
#
# rotate = TRUE, ggtheme = theme_pander(),
# title = "Distribution of connections") +
# scale_y_continuous(limits = c(0,15), breaks=c(0,5, 10, 15)) +
# scale_fill_discrete(breaks=rev(df.c2$NetworkName), ) +
# scale_fill_discrete(guide = guide_legend(reverse=TRUE)) +
# theme(legend.position = "right",
# axis.text.y = element_blank(),
# axis.title.y = element_blank(),
# plot.title = element_text(size = 20),
# axis.text = element_text(size = 20),
# legend.title = element_text(size = 20),
# legend.text = element_text(size = 20))
if (SAVE_PLOT) {
png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure10.png", bg = "transparent",
width =3000, height = 1000,
units = "px", res = 150)}
ggarrange(p7, NULL, NULL, NULL, p8,
#labels = c("A", "B", "C"),
ncol = 5, nrow = 1, align = "h",
common.legend = TRUE, legend = "bottom", widths = c(1,.8,.8,1))
Chi-sq Test
chisq.test(df.c1$count, df.c2$count, simulate.p.value = T, p = noi_stats$N)
##
## Pearson's Chi-squared test with simulated p-value (based on 2000
## replicates)
##
## data: df.c1$count and df.c2$count
## X-squared = 98, df = NA, p-value = 0.2994